You can find the slides (and additional materials) here: https://github.com/jstbcs/ws-bayesian-stats-r.

Julia Haaf

Assistant Professor at the Psychological Methods Department, University of Amsterdam

email: j.m.haaf@uva.nl

Job alert: Two temporary assistant professor positions at the Psychological Methods Unit (https://www.academictransfer.com/nl/297540/assistant-professor-of-psychological-methods/).

1. Introduction to Bayesian inference
2. The BayesFactor package
3. \(t\)-Test
4. ANOVA
5. Linear model setup

• Effect of media stories on attitudes towards Syrian refugees.
• IV: Sympathetic media story about refugees or control story about food trends.
• DV: Participants’ favorability ratings towards refugees.
• Also recorded: Political affiliation of the participants (conservative or liberal).

Example from Rouder, Haaf, & Aust (2018)

• Conduct a \(t\)-test to assess the effect of media stories on attitudes towards Syrian refugees.
• Conduct a 2x2 ANOVA to assess the effects of both media stories and political affiliation as well as their interaction.
• Not possible to provide evidence against these effects using classical statistics.
• Neither necessarily match the research questions associated with the study.
• Baysian model comparison allows to better match inference and theoretical questions.
• But you can also do the regular stuff.

1. The plight of refugees generates an empathic response which results in an increase in favorability for conservative and liberals alike.
2. Confirmation bias (Tversky & Kahneman, 1974): People ignore media stories incompatible with their world view. Therefore, only liberals are affected.
3. Belief polarization (Cook & Lewandowsky, 2016): Liberals see refugees more favorably after viewing positive stories while conservatives might actually see them more negatively.

• We need models

• We need models
• \({\cal M}_1:\) The effect of media stories on attitudes is positive.

• We need models
• \({\cal M}_1:\) The effect of media stories on attitudes is positive.
• \({\cal M}_0:\) There is no effect of media stories on attitudes.

• The prior distribution indicates the plausibility of effects before the data have been collected.
• Because we make this fine-grained specification, we can also make predictions before collecting any data!

• In Bayesian statistics, evidence for models reflects how well they predict data.
• Prior is on true values.
• Predictions are on data.

• In Bayesian statistics, evidence for models reflects how well they predict data.

• In Bayesian statistics, evidence for models reflects how well they predict data.

• In Bayesian statistics, evidence for models reflects how well they predict data.

• In Bayesian statistics, evidence for models reflects how well they predict data.

• Based on Rouder, Speckman, Sun, Morey, & Iverson (2009) (\(t\)-test), Rouder et al. (2012) (ANOVA), and Rouder & Morey (2012) (regression), total of > 4000 citations.
• Same setup is implemented in JASP and SPSS.
• Provides a convenient interface for commonly used models.
• Bayesian versions of typical frequentist tests using Bayes factor.

• General linear models (including linear mixed effects models): generalTestBF, lmBF
• Linear regression: regressionBF, lmBF
• Linear correlation: correlationBF
• \(t\)-tests: ttestBF
• Meta-analytic \(t\)-tests: meta.ttestBF
• ANOVA: anovaBF, lmBF
• Contingency tables: contingencyTableBF
• Single proportions: proportionBF

• General linear models (including linear mixed effects models): generalTestBF, lmBF
• Linear regression: regressionBF, lmBF
• Linear correlation: correlationBF
• \(t\)-tests: ttestBF
• Meta-analytic \(t\)-tests: meta.ttestBF
• ANOVA: anovaBF, lmBF
• Contingency tables: contingencyTableBF
• Single proportions: proportionBF

Can be used to test our position 1: The media story increases favorability ratings of refugees.

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study)

If you want to understand formulas in R better try this tutorial (for general setup) and this tutorial (for mixed models).

• \({\cal M}_1:\) The effect of media stories on attitudes is positive.
• \({\cal M}_0:\) There is no effect of media stories on attitudes.

• Attitudes might be normally distributed: \(Y_{ij} \sim \mbox{Normal}(\mu + x_j \theta, \sigma^2)\).
• For model \({\cal M}_0\) \(\theta = 0\).
• For model \({\cal M}_1\) we need a prior distribution on the effect, \(\theta\)
• Prior on \(\sigma^2\) and \(\mu\)?

1. It is necessary
2. You gotta know something

1. It is necessary
2. You gotta know something
3. Informed priors \(\rightarrow\) informative tests

Obtain strong evidence by increasing model discriminability.

Stefan, Gronau, Schönbrodt, & Wagenmakers (2019)

1. It is necessary
2. You gotta know something
3. Informed priors \(\rightarrow\) informative tests
4. Informed priors ensure identifiability

• Complex models often suffer from identifiability issues.

• Informed priors introduce probabilistic constraints.

• BayesFactor package in R and JASP use default priors for the t-test called JZS prior.
• Prior structure is consistent.
• Joint prior on \(\mu\) and \(\sigma^2\) is a Jeffreys prior.
• User-defined informed prior settings are on the scale of the effect size:
• \(\delta = \theta / \sigma\) (think Cohen’s \(d\)).

• BayesFactor package in R and JASP use default priors for the t-test called JZS prior.
• Prior structure is consistent.
• User-defined prior settings are on the scale of the effect size.

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study
                     , rscale = 1 / sqrt(2))

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study
                     , rscale = 1 / sqrt(2))

par(mar = c(3,3,.5,.5), mgp = c(2,.7,0))
x <- seq(-3, 3, .01)
y <- dcauchy(x, 0, 1/sqrt(2))
plot(x, y, type = "l", lwd = 2
     , ylab = "Density", xlab = "Effect Size"
     , col = "darkblue", ylim = c(0, .5))
lines(x, dnorm(x), col = adjustcolor(1, .5), lty = 2, lwd = 2)
abline(v = c(-1, 0, 1), lwd = 0.5, col = adjustcolor(1, alpha.f = 0.3))

• Prior \(\rightarrow\) we can make predictions on data.
• This is a bit tricky with the JZS prior though…

sim.pred.fun <- function(scale){
  delta <- rcauchy(1, location = 0, scale = scale)
  mu <- runif(1, -1000, 1000)
  s2 <- MCMCpack::rinvgamma(1, 0.5, 0.5)
  
  I <- 100
  cond <- rep(c(-1/sqrt(2), 1/sqrt(2)), each = I)
  y <- rnorm(2 * I, mu + cond * delta * sqrt(s2), sqrt(s2))
  
  return(diff(tapply(y, cond, mean)))
}

res <- replicate(100000, sim.pred.fun(1/sqrt(2)))
round(quantile(res), 2)

##          0%         25%         50%         75%        100% 
## -1179562.46       -1.74       -0.01        1.70   155159.42

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study
                     , rscale = 1 / 4)

1. Plot the default prior on effect size.
2. Adjust the scale until it matches your expectations of the outcome of the study.
3. Conduct a \(t\)-test using your setting.
4. Do you think any two people in the workshop get the same result?

Put your setting in the zoom chat!

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study
                     , rscale = 1 / 2)

## Bayes factor analysis
## --------------
## [1] Alt., r=0.5 : 2.142883 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS

The plight of refugees generates an empathic response which results in an increase in favorability for conservative and liberals alike.

• Does the \(t\)-test as we used it address this statement?

est.attitude <- BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study, rscale = 1 / 2
                     , posterior = T, iterations = 20000)

est.attitude <- BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study, rscale = 1 / 2
                     , posterior = T, iterations = 20000)

plot(est.attitude[,"beta (Neutral - Refugee plight)"])

The plight of refugees generates an empathic response which results in an increase in favorability for conservative and liberals alike.

• Does the \(t\)-test as we used it address this statement?

• We need an ordinal constraint.
• Ordinal constraint: Specification of the direction of an effect.

• One-sided \(t\)-tests are easy with BayesFactor.

• One-sided \(t\)-tests are easy with BayesFactor.

BayesFactor::ttestBF(formula = attitude ~ story
                     , data = media_study
                     , rscale = 1 / 2
                     , nullInterval = c(-Inf, 0))

## Bayes factor analysis
## --------------
## [1] Alt., r=0.5 -Inf<d<0    : 4.196947   ±0%
## [2] Alt., r=0.5 !(-Inf<d<0) : 0.08881904 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS

1. New data set, this time real data.
2. Run a Bayesian t-test, try out ordinal constraints and estimation, maybe even prior prediction.
3. Ask questions while I am here!

SourceURL <- "https://raw.githubusercontent.com/PerceptionCognitionLab/data0/master
/contexteffects/FlankerStroopSimon/cleaning.R"
devtools::source_url(SourceURL)

library(tidyverse)

stroop.agg <- stroop %>%
  group_by(ID, congruency) %>%
  summarize(mrt = mean(RT)) %>%
  spread(congruency, mrt)

• Based on Rouder et al. (2012).
• Between-subject, within-subject and mixed designs.
• Random effects.
• Perhaps the most thought out functionality in BayesFactor.

Figure from Rouder, Engelhardt, McCabe, & Morey (2016)

• Extension to \(t\)-test.
• Effects are parameterized as \(\delta = \theta / \sigma\) (effect size).
• Coding is orthonormal (depends on the number of levels per factor).
• Most simple case is two levels: \(x_1 = -1/\sqrt{2}, x_2 = 1/\sqrt{2}\).

media_study$story <- factor(media_study$story)
media_study$political_affiliation <- factor(media_study$political_affiliation)

BayesFactor::anovaBF(formula = attitude ~ story * political_affiliation
                     , data = media_study)

rscaleEffects = A named vector of prior settings for individual factors, overriding rscaleFixed and rscaleRandom. Values are scales, names are factor names.

BayesFactor::anovaBF(formula = attitude ~ story * political_affiliation
                     , data = media_study
                     , rscaleEffects = c("story" = 1/2
                                         , "political_affiliation" = 1/2
                                         , "story:political_affiliation" = 1/3))

BayesFactor::anovaBF(formula = attitude ~ story * political_affiliation
                     , data = media_study
                     , rscaleEffects = c("story" = 1/2
                                         , "political_affiliation" = 1/2
                                         , "story:political_affiliation" = 1/3))

## Bayes factor analysis
## --------------
## [1] story                                                       : 1.835452 ±0%
## [2] political_affiliation                                       : 191.5876 ±0%
## [3] story + political_affiliation                               : 468.3391 ±0.73%
## [4] story + political_affiliation + story:political_affiliation : 1156.833 ±7.65%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

• How does your prior influence the results of the analysis?
• Redo the analysis for a range of priors.
• If the results are (relatively) stable then we may trust them more.

• How does your prior influence the results of the analysis?
• Redo the analysis for a reasonable range of priors.
• If the results are (relatively) stable then we may trust them more.

Explore the range of settings using all combinations.

• In a multiverse analysis you draw your conclusions based on the entire space of possible parameter settings.
• The sensitivity analysis provides information about the robustness of your analysis and conclusions, it is not meant to replace the main analysis.
• Remember, you set the initial prior with care.

1. Run the ANOVA analysis on the dat_media.csv.
2. Conduct a sensitivity analysis for all three prior scales using reasonable ranges of settings.
3. Is the variability in Bayes factor reasonable?

geturl <-"https://raw.githubusercontent.com/jstbcs/ws-bayesian-stats-r/
main/data/dat_media.csv"
media_study <- read.csv(geturl)

• Not as straight-forward as with the \(t\)-test.
• Not as impossible as with frequentist ANOVA.

1. The plight of refugees generates an empathic response which results in an increase in favorability for conservative and liberals alike.
2. Confirmation bias (Tversky & Kahneman, 1974): People ignore media stories incompatible with their world view. Therefore, only liberals are affected.
3. Belief polarization (Cook & Lewandowsky, 2016): Liberals see refugees more favorably after viewing positive stories while conservatives might actually see them more negatively.

1. The plight of refugees generates an empathic response which results in an increase in favorability for conservative and liberals alike.
2. Confirmation bias (Tversky & Kahneman, 1974): People ignore media stories incompatible with their world view. Therefore, only liberals are affected.
3. Belief polarization (Cook & Lewandowsky, 2016): Liberals see refugees more favorably after viewing positive stories while conservatives might actually see them more negatively.

1. Belief polarization (Cook & Lewandowsky, 2016): Liberals see refugees more favorably after viewing positive stories while conservatives might actually see them more negatively.

\({\cal M}_2:\) The effect of the refugee media story is positive for liberals and negative for conservatives.

1. Belief polarization (Cook & Lewandowsky, 2016): Liberals see refugees more favorably after viewing positive stories while conservatives might actually see them more negatively.

\({\cal M}_2:\) The effect of the refugee media story is positive for liberals and negative for conservatives.

Encompassing approach (Klugkist, Laudy, & Hoijtink, 2005):

1. Calculate/estimate the prior probability \(P(\mbox{constraint})\) for the ordinal constraint.
2. Estimate the posterior probability \(P(\mbox{constraint} | Y)\) from posterior samples of the model.
3. The Bayes factor between the unconstrained and the ordinal-constrained model can be estimated as \[BF = \frac{P(\mbox{constraint})}{P(\mbox{constraint} | Y)}\]

• Two constraints, symmetric priors: \(0.5 \times 0.5 = 0.25\).
• Can also be estimated: \[ \begin{align} Y_{RL} &= \mu + \theta_1 + \theta_2 + \theta_3\\ Y_{CL} &= \mu - \theta_1 + \theta_2 - \theta_3\\ Y_{RC} &= \mu + \theta_1 - \theta_2 - \theta_3\\ Y_{CC} &= \mu - \theta_1 - \theta_2 + \theta_3\\ \end{align} \]
• Constraints refer to \(Y_{RL} - Y_{CL}\) and \(Y_{RC} - Y_{CC}\), respectively.
• \[ \begin{align} Y_{RL} - Y_{CL} &= 2\theta_1 + 2\theta_3\\ Y_{RC} - Y_{CC} &= 2\theta_1 - 2\theta_3\\ \end{align} \]

\[ \begin{align} Y_{RL} - Y_{CL} &= 2\theta_1 + 2\theta_3\\ Y_{RC} - Y_{CC} &= 2\theta_1 - 2\theta_3\\ \end{align} \]

M <- 100000
theta1 <- rcauchy(M)
theta3 <- rcauchy(M)

mean(theta1 + theta3 > 0 & theta1 - theta3 < 0)

## [1] 0.24943

tmp <- BayesFactor::anovaBF(formula = attitude ~ story * political_affiliation
                     , data = media_study
                     , rscaleEffects = c("story" = 1/2
                                         , "political_affiliation" = 1/2
                                         , "story:political_affiliation" = 1/3))

est <- BayesFactor::posterior(model = tmp
                              , iterations = 20000
                              , index = 4)

Reconstructing cell means and effects from the estimates.

Y_cl <- est[, "story-Neutral"] + 
  est[, "story:political_affiliation-Neutral.&.Liberals"]

Y_rl <- est[, "story-Refugee plight"] + 
  est[, "story:political_affiliation-Refugee plight.&.Liberals"]

Y_cc <- est[, "story-Neutral"] + 
  est[, "story:political_affiliation-Neutral.&.Conservatives"]

Y_rc <- est[, "story-Refugee plight"] + 
  est[, "story:political_affiliation-Refugee plight.&.Conservatives"]

(postprob <- mean(Y_rl - Y_cl > 0 & Y_rc - Y_cc < 0))

## [1] 0.36585

Bayes factor:

(BF_u3 <- 0.25 / postprob)

## [1] 0.6833402

bftab <- tmp@bayesFactor[, 1:2]
bftab[, 1] <- exp(bftab[, 1])
bftab[5, 1] <- bftab[4, 1] / BF_u3
rownames(bftab)[5] <- "Position 3"
knitr::kable(bftab)

• Go through the code and make sure you understand what is going on.
• Say we wanted to test whether the effect of media story is positive both for liberals and conservatives. What would change and what would stay the same?
• Implement the new constraint and report the Bayes factor for the new position 4 vs. the intercept-only model.

• General linear models (including linear mixed effects models): generalTestBF, lmBF
• Linear regression: regressionBF, lmBF
• Linear correlation: correlationBF
• \(t\)-tests: ttestBF
• Meta-analytic \(t\)-tests: meta.ttestBF
• ANOVA: anovaBF, lmBF
• Contingency tables: contingencyTableBF
• Single proportions: proportionBF

This function computes Bayes factors corresponding to restrictions on a full model.

data(puzzles)
BayesFactor::generalTestBF(RT ~ shape * color + ID
                           , data = puzzles
                           , whichRandom = "ID")

## Bayes factor analysis
## --------------
## [1] shape                            : 0.6114754 ±0.01%
## [2] color                            : 0.6114754 ±0.01%
## [3] shape + color                    : 0.3865341 ±0.82%
## [4] shape + color + shape:color      : 0.1328641 ±1.55%
## [5] ID                               : 111516.6  ±0%
## [6] shape + ID                       : 316107.6  ±0.98%
## [7] color + ID                       : 313479.8  ±1.41%
## [8] shape + color + ID               : 1309516   ±1.67%
## [9] shape + color + shape:color + ID : 512538.7  ±1.99%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

This function computes Bayes factors corresponding to restrictions on a full model.

This function computes Bayes factors, or samples from the posterior, of specific linear models (either ANOVA or regression).

BayesFactor::lmBF(RT ~ shape + color + shape:color + ID
     , data = puzzles
     , whichRandom = "ID")

## Bayes factor analysis
## --------------
## [1] shape + color + shape:color + ID : 488891.2 ±3%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

Rouder & Morey (2012)

• Now it is time to play.
• Repeat the analysis for dat_media using generalTestBF and lmBF.
• Alternatively, use your own data or continue with bugs.csv.

• BayesFactor can do a lot of conventional tests, but in Bayesian.
• Model comparison instead of parameter testing.
• For a full Bayesian analysis there are many tools in addition to Bayes factor:
• Prior predition,
• Estimation,
• Sensitivity analysis,
• Data plotting(!), …